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SUMMARY 


An algorithm for laminar and turbulent viscous compressible 
two'-dimensional flows is presented. For the application of precise boundary 
conditions over an arbitrary body surface, a body-fitted coordinate system is 
used in the physical plane. A thin»-layer approximation of the Navier-Stokes 
equations is introduced to keep the viscous terms relatively simple. The 
flow field computation is performed in the transformed plane. A factorized, 
implicit scheme is used to facilitate the computation. Sample calculations, 
for Couette flow, developing pipe flow, isolated airfoil, 2-D compressor 
cascade flow, and segmental compressor blade design are presented. To a 
certain extent, the effective use of the direct solver depends on the user's 
skill in setting up the gridwork, the time step size and the choice of the 
artificial viscosity. The design feature of the algorithm, an iterative 
scheme to correct an assumed geometry for a specified surface pressure 
distribution, works well for subsonic flows. A more elaborate correction 
scheme is required in treating transonic flows where local shock waves may be 
involved. 
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VISCOUS COMPRESSIBLE FLOW DIRECT AND INVERSE COMPUTATION WITH ILLUSTRATIONS 


1 . Introduction 
1 . 1 Background 

In the late 1960's, we extended Stanitz's (1) inverse solution method 
for planar potential flows to axisymmetric flows. Curved-wall diffusers were 
designed by using the extended method, and they are referred to as Griffith 
diffusers (2,3). In Griffith diffusers, the measured wall pressure 
distributions correlate very well with the prescribed distributions. In the 
early 1980' s, we extended the design procedure to include shear flows to 
design curved-wall diffusers for short ejectors (4). Even though record high 
thrust augmentation ratios were observed, the diffuser pressure distributions 
no longer correlated as well as they did in the sixties. The major 
difference is that in the Griffith diffusers, there is always a distinct 
potential core flow while in the ejector diffusers, the entire flow field is 
often dominated by the viscous forces. It is a natural progression that we 
began to take on the development of an inverse solution method for flows 
without neglecting viscosity. The momentum equation considered here-in is 
the Navier-Stokes equation with the thin-layer approximation. The energy 
equation is also included as one of the governing equations to account for 
the compressibility of the fluid. In the Reynolds stress terms, both 
Baldwin-Lomax (5) and Cebeci's (6) turbulence models were considered. The 
Baldwin-Lomax model was used in the cascade computations and Cebeci's model 
was used in the pipe flow computations. The NASA Lewis Research Center 
initially was interested m the possibility of applying this inverse 
procedure for gas turbine blade redesign. As the algorithm and computer code 
were being developed, its interest was shifted to the applicability of the 
code to compressor blade design. Therefore, a compressor blade of known 
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surface pressure distribution, reported by Schmidt et al (7), was used to 
verify the algorithm and the code. 

1 . 2 Outline of the Computational Procedure 

The computation begins with an assumed geometry. The flow field 
computation, including pressure distributions over the boundaries, is 
performed in a transformed plane. The transf ormation is from a body-fitted 
coordinate system. Thompson’s (8) general coordinate transformation is used 
in conjunction with Sorenson's (9) method to provide the orthogonality of the 
grid work at the solid boundaries. An implicit scheme for solving the 
compressible Navier-Stokes equation was developed, by the utilization of 
Beam and Warming's (10) scheme. The direct solver used in this report 
differs from Steger's (11) in that the present procedure allows solutions of 
axisymmetnc flows and planar flows while Steger’s is for planar flows only. 
The method of obtaining the desired geometry of either an axisymmetric flow 
passage or a blade cascade is achieved by the Secant method (12), and the 
virtual velocity method (13). These are progressive, iterative procedures. 
Thompkins and Tong (13) used an iterative method in correcting geometries, 
for inviscid flows. Modifications and interfacing of the above 
techniques were required to yield this comparatively complex procedure. In 
this report, we consider planar flows as sub- cases of axisymmetnc flows. In 
the sub- cases, appropriate terms are deleted from the axisymmetric flow 
equations. Flags are used in our computing code to direct the computation to 
either axisymmetric flow or planar flow. The listing of the computer code is 
given in Ntone's dissertation (1H) and will be deposited at Cosmic. Requests 
for copies of the code can be directed to Cosmic, Computer Software 
Management and Information Center, Suite 112, Barrow Hall, Athens, GA 30602. 
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2. Analysis 


In this chapter, the form of the Navier-Stokes equations selected for 
computation will be presented. The numerical algorithm selected for their 
solution will then be derived. A brief discussion on boundary conditions 
will be given here and additional details will be provided in the chapter of 
sample calculations. A brief presentation of Cebeci's and Baldwin-Lomax 
turbulence models for Reynolds stress calculation will be given. For 
improving computational stability at high Reynolds numbers, an artificial 
viscosity was introduced. It will also be discussed in this chapter. 

2. 1 Governing Equations 

For a viscous, compressible axisymmetric flow, the Navier-Stokes 
equations are written in the vectorial form similar to Peyret and Viviand 
(15). as 
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The above equations are written in the physical plane with cylindrical 

coordinates (x, r, t) for axisymmetric flows, where x denotes the axial 

coordinate, r the radial coordinate, and t the time. If the multiplication 

1 3 

factor r is set equal to one, the — terms go to zero, and the terms are 

g 

replaced by -5— then the sygtem of equations becomes useful for planar flows, 
dy 

In the above, p is the density of the fluid, u and v the velocity components 
in x and r direction, respectively, e the sum of thermal energy C V T and the 
kinetic energy (u 2 + v 2 ), p the pressure, Y the ratio of specific heats and 
p the viscosity. The equations are non-dimensionalized using reference 
values of length, velocity, density and viscosity. The Prandtl number Pr is 
defined by 


Pr = 


1J ref C p 

k ref 


where C p is the specific heat at constant pressure, and kpef * s a reference 
thermal conductivity. The Reynolds number, Re, is defined by 

Re ■ (-^) , 

v u ref 
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where L, is the reference length. Pressure, p, is obtained from the equation 


of state 

p = (Y-1 ) [e - 1/2 p ( u 2 + v 2 )] 


For laminar flow, the dependence of y on temperature can be accounted for by 
using Sutherland's relation: 


(T 

T 


,3/2 


ref 


C 2 + T ref 
C 2 + T 


where T denotes temperature, and 
C 2 = 198.6 °R 

For turbulent flows, it is customary to express each dependent variable as 
the sum of a mean and a fluctuating quantity, such as <j> = <(> + <}>' , and then 
time-average the Navier-Stokes equations. For compressible flows, this 
procedure leads to the presence of second and third order moments of the 
fluctuating variables due to the density fluctuations. To avoid this, the 
concept of mass averaging (16-18) is introduced. For example, if <j> 
represents the instantaneous value of a dependent variable, then the 
following decomposition is used: 

<(> = <j! + <j>" 

where $ is a mass averaged quantity, defined as $ = p4>/p, <f>" is a fluctua- 
tion, and the <j> denotes the time averaged quantity. The <t>" is then 
related to <j>', which is the customary fluctuation quantity, by 



P 


For Navier-Stokes equations, the mass averaged decomposition is applied to u 
and v while the customary time averaging is applied to p, e and p. The time 
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averaged turbulent Navier-Stokes equations are essentially the same as those 
for laminar flows, except that one must add to each laminar shear stress and 


heat flux term its corresponding turbulent contribution resulting from 

fluctuations. Thus, for turbulent flows, a turbulent viscosity pt 13 added 

to the molecular viscosity p, and the coefficient ^ is replaced by 

Y (“p + pp ) * where Pr fc is a turbulent Prandtl number. However, proper 
t 

closure of the equations is not achieved until a method of calculating the 
values of p^ and Pr^ is introduced. This subject will be dealt with in the 
section on "Turbulence Model." 

In many problems, one wishes to have boundary conditions satisfied 
exactly on an arbitrarily shaped body surface. Therefore, the need for a 
coordinate transformation from the cartesian coordinate system (x, r, t) to a 
more general curvilinear system ( 5 , n, t') arises. With the coordinate 
transformation 

5 = 5(x, r, t) 
n *> n(x, r, t) 

V - t , 

the axi3ymmetric Navier-Stokes equations can be rewritten as: 
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G = n. q + n F/J + n G/J 
t n x p 


G = n f,/j + n G./j 
1 x 1 r 1 


H = H/J 


J = £ n - £ n 
^x r r x 


1 /( Vn' Vc’ 


If one intpoduces the velocities, V, along n = constant line and, U, along £ 
= constant line, then F and G become: 

~ P U 

_ puU + P 

F = — 

J pvU + P 

(e + p) U - C fc p 


PV 

puV + n x p 
pW + n r p 
( e + p)V - n fc p 
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U £ «t + u ^x + v ^p 



t 


V s \ + un x + vn p • 

In the above tpansfopmed equations, the presence of the matrices £ x , Cp> 
n x , n r and the Jacobian, J, implies that the coordinate transf ormation is 
known. The details of how the transformation was obtained are given in 
Appendix A. For a stationary boundary, ^ and ht zero. For a moving 
boundary E,t and nt roust be specified. 

The viscous part, (RHS of Equation (2)), contains a large number of 
terms. These items can be greatly reduced if one uses the "thin layer" 
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assumption, according to which 5 -deri vati ves in the viscous part are 
neglected as small compared to the n~ derivative. Such an approximation is 
different from that of the usual boundary layer equation assumption, in the 
sense that it allows a pressure gradient across the streamlines even inside 
the boundary layer. Degani and Steger (19) recently presented a comparison 
of calculations using the thin layer assumption and the full Navier-Stokes 
equations to validate the approximation. When the thin layer approximation 
is incorporated into the governing equations, Equation (2) becomes: 


3q + II + 12 = _1_ 30^ + T7 
35 3n Re 3n 
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If we let the multiplication factor r be one in Equation (3), 1/r terms 
be zero, then H will vanish, and Equation (4) the governing equation for 
planar flows will be 
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The axisymmetric equations can also be written in a form very similar to 
the planar flow governing equations, and is shown as follows: 
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In this report examples are presented to show the solutions for 
Equations (3) and (4). Equation (5) is suited for more complicated 
geometries than equation 3. However, the computing time will be longer when 
Equation (5) is used. 

2.2 Numerical Schemes 

In this section the derivation of the finite difference equations used 
in the solution of Equation (3) Is presented. Following a similar procedure, 
the corresponding finite difference equations for Equations (4) and (5) were 
also derived but are not included in this report. 

The algorithm selected here is the implicit factored scheme of Beam and 
Warming. As a starting point, the time differencing formula is written in 
Pade's form as, 
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where A and V are forward and backward difference operators, and superscript 
"n" represents the nth time step. This is demonstrated by: 
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and t' = n(At') with positive integer of n. Recognizing Aq n ~^ = Vq n , 
Equation (6) can be rewritten as 
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-n , sat: 4( ^J , it: |a2 . 
ot ' 


1+C 


1+C 3t' 1+c 


+ o (e - e - j) (At* 2 ) + o (At' 3 ) 


(7) 


Using Equation (3) we can show that Equation (7) takes the following form: 
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2.2.1 Linearization 

Equation (8) is non-linear because AF n , AG n , AG n , and AH n 
contain terms depending on q n+1 . However, those terms can be linearized, 
as shown below for AF n : 
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From the definitions of F and G, it can be shown that: 
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G(q)/J = G(q/J) = G(q), 
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One may also show from matrix multiplication, that 
F(q) = Aq and G(q) = Bq. 

The Jacobian matrices, A and B, are given in Appendix B. Expanding F(q) 
in a Taylor series yields: 
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A similar expression for AG n can also be obtained, if £ is replaced by n. 
Equation (8) can be now expressed as follows: 
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In the present research, the three-point backward implicit scheme i 3 
used, where 0 = 1 and £ = V 2 . In the direct solver, the boundaries are fixed. 
Therefore, the metrics do not change with respect to time. It can then be 
shown that 

F n = F n and 5" = G n , 

and F n and G n no longer depend on (q) n+1 . One can approximate AG^ in 
the following form: 

.pH . r .7 tt \ n . -n 
AG 1 = (— M + M 1 ) Aq , 

and, similarly, 

AH" - (D + E) n Aq" . 


These are linear because M, M x , D and E are expressed at time n, and 
are therefore known, just as in the case of AF n and AG n . The matrices 
M lf M, D and E are listed in Appendix B. Then the linearized algorithm 
appears as f oil ows : 
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2.2.2 Factorization 

A major gain in computational efficiency may result if Equation (9) is 
written in the following factored form: 
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= R.H.S. of equation (10). 
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and 


— n+1 — n A — n 

q = q + Aq 

Equation (11), the factored equation, differs from Equation (10) only by 
terms which do not affect the stated accuracy. Using Equation (11), the 
two-dimensional operator of Equation (10) was reduced to two one-dimensional 
operators . 

For spatial discretization, the central difference scheme is used, and 

the terms (-I 7 O and are represented by the following approximations: 

H i,j 3n 3n lfj 


3«D *i+ 1 .J - Vl.j 0 ,.-2. 

uh.i 2Ll + ° (A5 > 


and 
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3n 6 ¥n i,j 2 (An) z_ 


+ o(An ), 


When A£ and An are taken to be 1, Equation (11) becomes: 
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2 .1 i 1-max • 

The Kjj represents the right hand side of Equation (10), and Lj--], Qj , N^+i 
are 4x4 matrices. The solution of Equation (12) for Aqfj requires the 
inversion of a block tridiagonal matrix. This procedure is shown in Appendix 
C. The solution of Aq's, of Equation (11), which is referred to as the 
q-sweep, is obtained from the following: 


— n ^ — n , .. — n *n 

j-i 4C| i,j-i * V q u * Vi 4q i,j*i ' 4 q i.j 


(13) 


where 


V> 


9At 1 
1 + C 




L 

Re i.j-1 


Re M1 i,j-1 " E i,j-1 ) 



6At' ,-n „n 

1+C i,j + 1 Re i,j + 1 


— Ml n 
Re i,j + 1 


E iJ + 1 ) 


Q 


J 


I 


0At * . 1 n 
1+C ( Re K l,j 
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and 

^ 1 J £ Jmax “ 1 • 

The matrices W, Y, R are given in Appendix B. For each time step, the 
solution of the previous time step is used to calculate the R.H.S. of 
Equation (10) and the L.H.S. of Equations (12) and (13). The solutions of 
Equations (12), (13) and (11a) yield the flow field at interior points of 
the new time step. At the boundaries, either the flow variables are 
specified, or they are calculated from the known values at the interior 
points. To initiate the calculation, an initial set of values must be 
assumed. The three-point backward scheme uses three time levels; the third 
level is calculated from the previous two levels. Therefore, one could not 
start the calculation with just an initial solution at n=0. The trapezoidal 
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rule Is employed for the first time iteration, and the three-point backward 
scheme is used thereafter. 

2.3. Boundary Conditions 

The importance of treating the boundary conditions correctly has long 
been pointed out by authors such as Moretti (20). The 1981 Symposium on 
Numerical Boundary Condition Procedures (21) was a recognition of the 
importance of this aspect in numerical computations. Essentially, boundary 
conditions affect both the well-posedness of initial- boundary value problems 
and the stability of their numerical solution (22). 

At in-flow either the velocity components and the density (or 
temperature) or the stagnation temperature and pressure and the flow angle 
were specified. At solid boundaries, both u and v are zero and either a 
constant surface temperature or an adiabatic wall condition is specified. 
Fluid properties which are not specified at boundaries are obtained from 
solutions at interior points by extrapolation. At the outflow boundary, 
where p is specified, p, u and v are obtained by linear extrapolation. 
Details on the extrapolation are provided in Appendix D. 

2.4 Turbulence Model 

As stated earlier, a complete closure of the governing equations 
require that a turbulence model be provided. In this research, two 
turbulence models, namely Cebeci's algebraic model and the variation of it 
introduced by Baldwin and Lomax (5), were examined. 

Cebeci's model (6) is a two-layer model in which the calculation of the 
eddy viscosity, y^, depends on whether one considers points in the inner 
region or the outer region. In the inner region, the expression for y^ is 
based on Prandtl's mixing length theory: 
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‘■Vinner ' 1 I ' 

where y is the normal distance from the wall and i the mixing length. They 
are obtained from the Van Driest expression 

i “ ^y [1-exp (-pp-) ] , 

with 

= 0.-4 

A + = 26v . 

Pw 

Here, v is the kinematic viscosity, and i w and p w are the shear stress and 
density at the wall, respectively. 

In the outer region, pt i3 calculated from: 

( y t^outer ~ w o°kleb p 

where 

00 

Uq - k 2 / (Ug-u)dy , 


and 


'kleb 


[U5.5(|) 6 ] _1 


In the above, u e is the velocity at the edge of the boundary layer, 6 is the 
boundary layer thickness, y is the Klebanoff intermittency factor, and k 2 = 
0 . 0168 . 

The inner region extends from the wall to the crossover point y c , where 
(pt)inner = Uo* For y > yc* values of (ut)outer 31,6 applicable. One of the 
difficulties in applying Cebeci's model lies in the accurate determination of 
6. Baldwin and Lomax (25) pointed out that significant errors would result 
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from the inaccurate evaluation of 6 and modified Cebeci's model to avoid this 
calculation. 

In the Baldwin and Lomax model, the eddy viscosity for the inner region 
is obtained from: 

(^t dinner = P ^ | u | * 

where 


w 


3u _ 3v 
3r 3x . 


For the outer region, the following formula is used 
(ut^outer =1< 2 *^cp P F wake ^kleb (y) • 

where 
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use the smaller of the two. 


F k le„ (y) ■ 0*5.5 (C kleb y/y max ) 6 ]'' 

Ckleb “ C-3 C W |< = 0.25 

F ma x the raaximun of F(y) defined as 


F(y) = | to | l , 


and y max is the y-location where F(y) maximum occurs, u^f is given 


by 


dif 


( / u z + v z ) 


max 


The success of this model relies on the existence of a well-defined peak 
for the function F(y) and not having to calculate the boundary layer 
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thickness, 6. In the current cascade computation, the Baldwin-Lomax model 
was used in the computer program. No difficulty was encountered. 

2.5 Artificial Viscosity 

The algorithm described in section 2.2 works well for low Reynolds 
number calculations, but often requires additional damping terms for 
computational stability at high Reynolds number. The form of the damping 
terms used here is the same as Steger's (11). When these damping terms are 
added. Equation (11) becomes: 
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(14) 


In Equation (14), a fourth order dissipation term has been added to the right 
hand side of the equation, while second order terms are added to both 
implicit factors. Operators such a3 75 and A5 are backward and forward 
operators. In a central difference approximation, 
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are used for second and fourth order derivatives, respectively. The 
coefficients, e, and e 2 , are dissipation coefficients. 

Artificial viscosity terms of this nature were studied by Desideri et 
al. (22). Their study was based on simplified forms of the governing 
equations and certain types of boundary conditions. Nevertheless, their 
limits for the ratio of e 2 to are used as a guideline for the stability 
behavior of the present damped algorithm. 

2.6 Inverse Design 

The inverse design consists of computing a wall geometry that satisfies 
a specified pressure distribution. One may initiate the inverse design with 
an educated guess-geometry and a reliable procedure for wall modification. 
Each newly selected wall geometry should yield a wall pressure closer to the 
target pressure than the previous ones. In this research, two methods have 
been examined. One is an adaptation of the so-called secant method, often 
used in the numerical solution of non-lmear equations. The other is based 
on the virtual wall velocity idea previously used by Thompkins and Tong (13). 
2.6.1 Secant Method 

The secant method is illustrated in the following sketch. In this 
figure, it is desired to find the zero of a function of <p ( x ) . 



Starting from two initial guesses ( x -j , 4> -j ) and on the curve, a 

straight line is drawn through those two points, intersecting the x-axis at 
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X3. Then another straight line is drawn through (x2i<J> 2) and ( x 3»<l ! 3)» 
intersecting the x-axis at x^. If this process is continued, the 
intersection of each new straight line with the x-axis should get closer and 
closer to the correct value of x for which $ vanishes. If the two initial 
guesses are sufficiently good and the function <f> sufficiently well-behaved, 
convergence of the method is quite rapid. 

For the extension of the method to inverse design problems, the 
following formula is used: 



(15) 


In Equation (15), all quantities are evaluated at the wall. Subscripts "1" 
and "2" denote the two guesses used to calculate the next wall geometry and 
subscript "d" the design pressure. At each wall correction, Equation (15) 
must be used at every wall point in the region of interest. Clearly, this 
means smooth wall geometries will be obtained only if wall pressures are 
smooth. Since this may not necessarily be so, a wall geometry smoothing 
method must be introduced. In this case, all calculated wall geometries were 
smoothed by a least squares method, using a low degree polynomial. 

In order to apply the above wall correction scheme to a design problem 
using the thin layer Navier-Stokes equations, it is possible to perform a 
complete direct calculation on each of the two geometries used to obtain the 
next one and repeat the process until a satisfactory geometry is obtained. 
However, it seems computationally more efficient to treat the whole problem 
as a single direct computation in which wall corrections will be made after a 
specified number of iterations (time steps) of the direct solver. One may 
start the calculations using one guess for the geometry. The direct 
calculations are stopped when the maximum absolute difference between the 


2k 



current wall pressure and the wall pressure at the previous time step (wall 
pressure residual) becomes smaller than a specified value Op. At this point, 
the second guess-geometry is used, and direct calculations proceed with this 
second geometry until the wall pressure residual falls below Op. Then using 
the previous two geometries and their corresponding pressures, a new or the 
third geometry is calculated using equation (15). This geometry and one of 
the two guess-geometries whichever has a pressure closer to the design value 
are used to obtain the 4th geometry. The process repeats with 3rd and 4th 
geometries as the first and second guess geometries. For accuracy, a small 
value of Op should be specified. However, since o p determines the number of 
direct solver iterations between two wall corrections, the o p value may be 
increased to speed up the computation. For a given o p , the number of direct 
solver iterations decreases as one approaches the target geometry. 

Therefore, to maintain accuracy, one may decrease o p as the target geometry 
is approached. 

The form of Equation (15) makes it suitable for the design of flow 
passages such as nozzles, diffusers, etc. where the inlet radius is specified 
and remains fixed throughout the calculations. However, at the exit, where 
pressure is usually specified, the term P 2 _ Pi vanishes, which makes the 
formula inapplicable. This inconsistency affects points in the neighborhood 
of the exit, which must be obtained by extrapolation. Here, the 
extrapolation was done so that the flow can be made nearly parallel at the 
exit by maintaining the exit wall radius equal to the wall radius of the 
adjacent point. This, of course, imposes a limitation on the kinds of exit 
flows to be considered, but there does not seem to be a general rule of 
extrapolation for the exit wall nodal points. 
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Due to the presence of ?2 ~ Pi in the denominator in Equation (15), it 
is clear that numerical problems will arise if the difference between 
calculated and target pressures becomes very small. Therefore, calculations 
must be stopped when that difference becomes smaller than a pre- determined 
value. Also, initial guesses must be selected such that their corresponding 
pressure distribution curves do not cross each other. 

2.6.2 Virtual Wall Velocity Method 

Referring to the flux vector, G, at the right hand side of Equation 
(10) and imposing the no-slip condition at the wall, it can be shown that: 



at all times for the direct calculations (for a fixed wall). The subscript, 
w, refers to the location at the wall. However, it is possible to imagine a 
problem where the wall would be moving towards a steady state position 
corresponding to the target (design) pressure. At that steady state 
position, G w would be: 



where superscript "d" denotes the design conditions. Assuming that the wall 

moves with a finite velocity towards its steady state position, u, v, V and 

nt are no longer zero at the wall. They can be calculated by requiring G w 
for the moving wall to be equal to G^, which leads to: 
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Equation ( 1 6 ) represents a system of four equations. The first equation 
is automatically satisfied when the steady state is reached. If 
geometric changes in time are disregarded, the last equation yields: 


n t = + l) V, 


and the second and third equations yield: 
n x 

U “ V 

n r 

One may then obtain v from the second equation as: 
2 


v =■ + ( 


2 2 
n + n 
x r 


1 „ 1/2 
iip a -p|) • 


Now from the equation for nt and the definition of V, it can be shown that 

V = ' f (un x + v V * 

By applying the rules of partial differentiation, one can show that: 

n t ’ "Vf ‘ W • 

Assuming that the wall only moves in the r-direction, then 
n t n r r t’ ’ 

and the magnitude of the wall correction is given by: 


| Ar | = At 
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( 17 ) 
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If the calculations are initiated with an assumed geometry, Equation (17) can 
be used to correct the wall at each time step. However, since this would 
require the calculation of a new grid at each time step, it is useful to 
introduce Op as described in 2.6.1 to set the number of time steps between 
two wall corrections. In practice, it may be necessary to multiply the right 
hand side of Equation (17) by a constant X. The value of X should be 
selected small enough to keep the calculation stable. A value of 0.05 made 
the computations very slow and a value of 0.2 was satisfactory. Also, the 
computed r's at the wall are smoothed out by adding to the newly computed r's 
the expression: 


i+1 ,1 


- 2 r. . 

i.l 


* Vi,i> 


where j=1 denotes the wall, and k a positive constant. In the calculations 
performed here, k = 0.2 was used. 
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3. Results and Discussion 


In this chapter, some of the cases used to test the algorithm of the 
numerical solution described in the analysis chapter are presented. Since 
the inverse design feature of the computational program relies entirely on 
the direct solution, it is important to assess the capabilities and 
weaknesses of the direct solver. Therefore, a great part of the chapter will 
be devoted to the results of the direct solver. Then illustrations will be 
provided for the inverse design procedure. These illustrative cases have 
either known analytical or numerical solutions or experimental data. All 
calculations were performed on an IBM 308 1 . 

3 . 1 Direct Solution Results 

3.1.1 Laminar Flows 

3. 1.1.1. Flow Formation in Couette Motion 

The first case is the two-dimensional flow between two parallel plates 
at low Reynolds number. One of the plates is held fixed and the other one is 
suddenly set into motion with a constant velocity, u 0 , at a Reynolds number 
of 6.2. Such a flow was considered by Beam and Warming (23), when they 
extended their algorithm, originally derived for hyperbolic systems in 
conservation law form, to the Navier-Stokes equations. However, their 
governing equations were the full Navier Stokes equations in physical 
coordinates, while a coordinate transformation and a thin-layer approximation 
are used here. 
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Figure 1 : Flow formation in Couette motion 

The geometry of the problem is shown in Figure 1. The no-slip condition is 
applied at the walls (j=1 , and j=jn) along with an adiabatic wall condition. 
At the ^-boundaries (i=l and i=i m ), periodic boundary conditions are imposed 
as follows: 

Ql * J = qi rrf 1 >J 

q2 -j = Ql m>J 

Essentially, the periodic boundary conditions are used here to maintain the 
non-dependence of the flow field on the x-direction (or ^-direction) . The 
solution of Equation (4) for this case is compared to the exact solution (24) 
in Figure 2. The selected grid was uniform and had 6 points in the 
^-direction and 11 points in the q-direction. As expected, the solution 
showed no variations in the x-direction, and only one x-station needed to be 
considered. Due to the use of a relatively large time step, only 40 
iterations were required for convergence. Because of a large time step, the 
computed results of 5 and 10 iterations, labeled as 5n and lOn in Figure 2, 
are in poor agreement with the exact transient solution, while the exact 
steady state solution and the computed result labeled 40n are in good 
agreement. This is an indication that the steady state solution is 
independent of the time step and of the initial conditions. Better 
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gure 2. 


Flow formation in Couette motion. 


^ C. U. Numerical, constant wall temp 



Figure 3. Centerline velocity for the entry region of a pipe 



agreement with the exact transient solution is obtained with smaller time 
steps, but more iterations are required. 

3. 1.1. 2. Entry region of a circular pipe at Rer) ~ Mo 

The second case is the entry region of a circular pipe at low Reynolds 
number, with a uniform velocity u 0 of 100 ft/sec prescribed at the inlet. 
Friedman et al. (25) used a numerical solution of the incompressible 
Navier-Stokes equations to show how the Reynolds number affects the flow. In 
Figure 3, the present solution of the centerline velocity is compared to 
theirs for a Reynolds number of 40. The present solution (Eq. 3) was 
obtained with a 21x15 grid and converged in 400 iterations from an impulsive 
start from uniform conditions. Both a constant temperature wall and an 
adiabatic wall were considered. It is seen that the constant temperature 
condition yields a better agreement with Friedman et al.'s solution. 

3. 1.1. 3. Entry region of a circular pipe at Ref) = 3000 

The third case is the laminar flow in the entrance of a pipe at Re<j = 
3000. The uniform velocity prescribed at the inlet was again 100 ft/ sec. 
Because the entry region at this Reynolds number is quite long (above 300 
radii), only part of it was considered in the present calculations. A 41x18 
grid was used, and the solution converged in 1200 iterations from an 
impulsive start at uniform conditions. The solution is shown in Figure 4, 
along with Nikuradse's experimental curves (26) and the analytical results of 
Campbell and Slattery (27). Also shown is the only point from Reshotko's 
experimental data, as presented in Reference 27, that falls within the range 
considered here. Figure 4 shows that the agreement of the present solution 
with the solution by Campbell and Slattery (27) is generally good, except at 
the centerline. Nikuradse's data are also presented for comparison. The 
validity of these data for comparison in this problem, however, has been 
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— Nikuradse (experimental) 

— Campbell and Slatery ( Analytical ) 



Figure 4. Entry region of a circular pipe. Laminar flow velocity distributions. 


— Schiller (analytical) 



Figure 5. Entry region of a circular pipe. Laminar flow pressure distribution. 
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questioned by Campbell and Slattery, Shapiro (see discussion in Reference 27) 
and Sparrow et al . (28). The data point from Reshotko's experiment seems in 
good agreement with the present solution near the centerline. Figure 5 shows 
a comparison of the computed non-dimensional present pressure distribution 
versus the non-dimensional axial distance divided by the Reynolds number with 
the analytical results of Schiller reported by Tietjens (26), are in very 
good agreement with experimental data. It should be noted that the present 
results were obtained using Equation (3) and a rather coarse grid with the 
spacing in the x-direction taken as one diameter. For pipe flow problems, 
the solution of Equation (3) at high Reynolds numbers seems to allow such 
large spacings. However, making Ax much larger than one diameter would 
result in a computationally unstable situation. 

3. 1.1. 4. Isolated airfoil 

As a test of the algorithm for more complicated geometries, the laminar 
flow over an isolated NACA 0012 airfoil was examined. The calculations were 
performed at a Reynolds number of 10,000 with respect to the chord length and 
a free stream Mach number of 0.2. The grid around the airfoil was obtained 
using Sorenson's code (9) which allows specification of a uniform rpspacing 
and near-orthogonality at the wall. For the computations, a 79x3*1 grid of 
the C-type, as shown in Figure 6, was used. The distance between the wall 
and the nearest £-line was approximately 0.0011, and an elliptic shape outer 
boundary was specified at 7 chords away from the airfoil. A zero-degree 
incidence angle (measured from the x-axis) was selected, which means symmetry 
of the flow field existed about the x-axis. For a geometry such as shown in 
Figure 6, the elliptic outer boundary is the inflow boundary, while the 
vertical line at the right is the outflow boundary. At the outflow boundary, 
the pressure was specified, while p, u and v were extrapolated. The 
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Figure 6. Body-fitted coordinates, NACA 0012 airfoil. 
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extrapolation was done by taking 3/3n as zero, where n is the direction 

normal to the outflow boundary. This involves the inversion of three scalar 

tridiagonal matrices to obtain p, u and v at the boundaries. At the inflow 

boundary, two methods of specifying the boundary conditions were considered. 

The first consisted in specifying p, u and v and extrapolating for P by 
3P 

taking -^ = 0. The second consisted of specifying the stagnation pressure 
P 0 , and temperature, T 0 , and the flow inlet angle, then the left running 
Riemann invariant, given by: 


where c is the speed of sound and oi the total velocity (u> = /u 2 +v 2 ) , was 
linearly extrapolated from the interior to the boundary (see Reference 29). 
For adiabatic flows, w at the inlet can be obtained from: 

(Y-i) r~ + (MY+ncpTo - 2 (y-i )(r ") 2 ) 1/2 , 

W ' (Y+1 ) 

where all quantities are evaluated at the inlet. Then u and v are obtained 
from w using trigonometric relations and the flow inlet angle. Also, knowing 
w and T 0 , the inlet static temperature and speed of sound are obtained. The 
inlet Mach number can then be calculated, and the inlet pressure is obtained 
from an isentropic relation between P 0 and the Mach number. This method is 
valid for inviscid flows and only approximately correct for viscous flows. 
Therefore, the inlet boundary should be taken far enough from the airfoil. 

For the isolated airfoil problem, both methods of treatment of the inflow 
boundary were found adequate. 

At the wall, the no-slip condition (u=v=0) was specified, along with an 
adiabatic wall condition. The wall pressure was obtained from the interior 
pressure, which requires the inversion of a scalar tridiagonal matrix. At 
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the common boundary in the wake region, all flow variables were extrapolated 
and the averages of the extrapolates from above and below were used. The 
extrapolation procedure, again, assumed no gradients in the normal direction. 

The results of this calculation are shown in Figure 7, where the 
computed pressure distribution is compared to Mehta's incompressible solution 
as presented in reference 11. The present solution, obtained after 500 
iterations from an impulsive start from free stream conditions, agrees well 
with Mehta's solution. 

3.1.2 Turbulent flows 

3. 1.2.1. Flow in a circular pipe 

In order to assess the validity of the turbulence models presented in 
Chapter 2, a turbulent flow in a pipe was considered. Calculations were 
performed in a pipe section 5 radii long. A logarithmic velocity profile was 
prescribed at the inlet to simulate fully-developed conditions, and the 
objective of the test was to verify that such a profile would be maintained 5 
radii downstream. This test was performed at a Reynolds number, of 
110,000, with an average velocity , u, of 100 ft/sec. Both the Cebeci 
and the Baldwin-Lomax models were examined, and the velocity profiles 
obtained at the pipe exit using both models are compared to the logarithmic 
profile in Figure 8. A better agreement with the logarithmic profile is 
obtained with Cebeci 's model. It was also observed from the computed data 
that the Baldwin-Lomax model tended to underpredict the wall shear. For that 
reason, Cebeci 's model was used in pipe- flow computations. 

As a more severe test for turbulent-flow calculations, the flow in the 

entrance region of a pipe was considered. The calculations were performed 

for a Reynolds number of 388,000 and an inlet velocity of 200 ft/ sec using a 

Ar 

31x34 grid. The grid spacing near the wall was — = 0.0005, with r 

r n o 
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cp=p-p»/M.ui 


Mehta (imcompressible, analytical) 



Figure 7. Pressure coefficient distribution, NACA 0012 airfoil, zero 
degree incidence, M=0.2 Re =10,000, laminar. 



Figure 8. Fully developed flow in a circular pipe at Re^llOjOOO. 
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being the pipe radius. This grid was stretched exponentially towards the 
centerline, in order to reduce computational time. In the x-direction, a 
uniform increment of one diameter was used. Due to the small grid near the 
wall, the time step for this problem had to be taken small enough to allow 
stable calculations from an impulsive start at uniform conditions. However, 
to accelerate the computations, the time step was increased by a factor of 20 
after 100 iterations. The velocity results after 600 iterations are shown in 
Figure 9, along with the experimental data of Barbin and Jones (29), obtained 
at the same Reynolds number. The agreement is generally good, although the 
numerical solution (of Equation 3) indicates a slightly larger acceleration 
near the centerline. In Figure 10, the computed pressure distribution is 
compared to the experimental one. The experiment was performed on a pipe 
length of 45 diameters, while the calculations were performed on a pipe 
length of 30 diameters. For comparison, the reference pressure was 
arbitrarily taken as the exit pressure of the computations. In this case, 
the numerical solution deviates only slightly from the experimental data. 

This could partially be attributed to the imperfections of the turbulence 
model, and partly attributed to the difference in pipe length between 
computational and experimental cases. 

3. 1.2. 2. Compressor cascade 

The last case of the direct solution is the flow in a compressor cascade. 
The cascade blade geometry used in the calculations was the result of an 
analytical and experimental investigation by Schmidt et al. (7). Sorenson's 
code (9) was used to obtain the finite difference grid for the cascade. The 
grid is shown in Figure 11. It is a 79x34 grid of the C type, in which the 
vertical line on the extreme right represents the outflow boundary, and the 
vertical line at the extreme left, the inflow boundary. The outflow and 


39 



n/n 


— Barbin and Jones, experimental (Re d = 388,000) 
0 ] C.U. Numerical (Re d = 388,000) 


o 



r/r 0 = 0 

r/r c = 0.499 

r/r a = 0.749 


r/r 0 = 0.936 


o.o 


10.0 

x/d 


20.0 


30.0 


Figure 9. Entry of a circular pipe. Turbulent flow velocity distributions. 



Figure 10. Entry region of a circular pipe. Turbulent flow pressure 
distributions . 
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inflow boundaries are connected by two periodic boundaries maintained at a 
vertical distance equal to the cascade spacing from each other. It should be 
noted that the periodic boundaries are joined to the inlet boundary by 
circular arcs. For a clean separation between inlet and periodic boundaries, 
it seems desirable to keep the radius of those circular arcs as small as 
possible. However, it was found that trying to use a very small radius could 
make the grid generation calculations unstable. 

The spacing maintained between the inner-boundary and the nearest £-line 
was approximately 0.00056. The blunt trailing edge of the airfoil was 
changed to a sharp trailing edge in order to allow a smooth stretching of the 
grid in the wake region (^-direction) , without the need for additional points. 
This modification was not expected to significantly affect the solution away 
from the trailing edge. 

The boundary conditions for the cascade problem were treated as 

described in the case of the isolated airfoil, except for the periodic 

boundaries. At the periodic boundaries, extrapolation from the interior 

(using - 5 — = 0) was used for all variables. In order to maintain periodicity, 
dn 

the average values of the extrapolators from above and below were imposed at 
both boundaries. 

The calculations were performed at a Reynolds number of 170,000, which 
is roughly a tenth of the Reynolds number of the experimental investigation 
of Reference 7. It was assumed that the solution at the lower Reynolds 
number would still represent the essential features of the solution at the 
higher Reynolds number. Attempting a high Reynolds number calculation would 
have required the use of a finer grid near the wall, and, therefore, a much 
larger number of grid points. This, combined with the smaller time steps 
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associated with finer grids, would have led to a much more expensive 
solution. 

Convergence for the cascade problem was' found to be quite slow, with 
large oscillations of the dependent variables in time. In particular, when 
p, u and v were specified at the inlet and the pressure extrapolated, large 
pressure fluctuations were observed at the inlet. Since those fluctuations 
did not seem to die out fast enough, this treatment of the inlet was 
abandoned in favor of the specification of the stagnation pressure and 
temperature and the flow angle (35.3°). Although oscillations were present 
in this inlet treatment as well, they were not as large, and they died out 
faster. However, in this case, the inlet Mach number cannot be specified but 
is determined from the specified inlet stagnation quantities and outflow 
static pressure. For given inlet stagnation quantities, the exact inlet Mach 
number will be obtained if the exact downstream pressure is specified. Since 
the exact values were not available at the required locations, the 
experimental inlet Mach number (0.75) could not be obtained without going 
through a trial- and- error procedure. 

The results of the calculations shown in Figure 12 are for inlet Mach 
numbers of 0.734 and 0.77 measured at the inlet midpoint. Referring to 
Figure 12, comparison of the two computed Mach number curves suggests that 
for a given upstream stagnation conditions, a small change in downstream 
static pressure may result in a small change in inlet Mach number which in 
turn may have a significant effect on the solution. Therefore, one may 
expect that the non-uniform inlet Mach numberswill affect the solution. 

Since a uniform inlet Mach number only occurs when the inlet boundary is 
taken far enough from the airfoil, it is reasonable to expect that doing so 
might improve the results. Figure 13 shows a linear interpolated result for 
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Figure 12. Surface Mach number distributions. 



A comparison of surface Mach number distributions between 
interpolated numerical results and experimental data. 


Figure 13. 



Figure 12 for inlet Mach number of 0.75. It is compared to the experimental 
data. The wall Mach numbers plotted in these figures were obtained from an 
lsentropic relation between the inlet stagnation pressure and the wall static 
pressure. In general, Figure 13 shows that a good agreement exists between 
the computed and the experimental values with an inlet Mach number of 0.75. 
However, in all cases the computed peak Mach number seems too high. 

Over the pressure side, which is the lower surface, the agreement 
between computed and experimental Mach number is very good for the first 60$ 
of the chord. A maximum deviation of 10$ occurs around 80$ chord location. 
Over the suction side, the upper surface, deviations occur in the 15 to 40$ 
chord region. The experimental peak seems to occur at 25$ chord while the 
computed peak occurs at 30$ chord. The maximum deviation, however, only 
amounts to 7$. 

The flow compression takes place rapidly from 30$ chord onward over the 
upper and lower surfaces. Some accleration occurs over the lower surface 
toward the trailing edge. This flow behavior is shown both computationally 
and experimentally. Because of the rapid compression, it is suspected that 
there might be a possibility of flow separation over the upper surface toward 
the trailing edge. The velocity field plot, shown as Figure I 1 *, indicates 
that there is a small flow reversal and then a reattachment over that region. 

Experience gained with the cascade calculations shows that care must be 
exercised m obtaining an adequate grid for the computations. The 
distribution of points along the airfoil must be such that the geometry is 
adequately represented. 

Convergence for cascade direct computations may require 3000 to 5000 
iterations from an impulsive start. Acceleration of the computations 
requires the use of large time steps. The non-dimensional time step At' 
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Figure 14. Velocity field between compression blades. Enlargement over 

the upper surface near the trailing edge showing flow reversal 
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considered here varied from 0.001 to 0.02. However, large time steps require 
more smoothing for stability. Since excessive smoothing may lead to less 
accuracy, a careful balance must be kept in the choice of the time step and 
the smoothing coefficient (referred as ei , in section 2.5). Also, subsonic 
computations require less smoothing than computations involving transonic 
effects. It should be noted that upwind differencing has not been 
incorporated in the code. Therefore, the code is not quite ready for 
transonic computations. This may explain some of the discrepancies observed 
in Figure 13 where the wall Mach number exceeds 1. On the IBM 3081 and for a 
79x3*1 grid, roughly 4.2 seconds of CPU time are required for each iteration. 
3.2 Inverse Solution Examples 

The illustrations presented for the direct solution were designed to 
show that good results can be obtained for a variety of cases. For the 
design problem, it is assumed that the direct solver is good, and the success 
of the inverse methods presented earlier is judged by their ability to 
produce a wall geometry that yields a pressure distribution close to a target 
distribution. 

3.2.1 Secant Method 

An axisymmetric nozzle with a wavy wall geometry as shown in Figure 15 
was selected for sample calculations. The corresponding pressure 
distribution, as calculated by the direct solver, is shown in Figure 16. 

This pressure distribution, along with the two initial geometry guesses shown 
in Figure 15, were input to the design calculations. Figures 15 and 16 show 
the convergence histories for the wall geometry and the wall pressure, 
respectively. It is seen from those figures that a good design can be 
obtained with a small number of wall corrections. In this case, the 
calculations were started with uniform flow conditions prescribed on the 





Design 

1st wail correction 
2nd wall correction 
3rd wall correction 
4 th wall correction 


Figure 15. Wall geometry at various interation levels. Secant method 
applied to a nozzle at Re^ =24. 



Figure 16. Wall pressure at various iteration levels. Secant method 
applied to a nozzle at re d =24. 
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first guess-geometry. The value of Op was kept constant throughout the 
calculations, but was selected small enough to allow accurate results. The 
449 direct solver iterations required for the problem are not excessive, 
given that the direct solution for a single geometry would require about 200 
iterations. 

3.2.2 Virtual Wall Velocity Method 
3. 2. 2.1. Flow nozzle 

The nature of the virtual wall velocity method restricts its use to 
problems in which the region to be corrected lies between two fixed points. 
For example, its use for a nozzle design problem may assume that the 
locations of the inlet and exit radii are fixed, and the interior part is to 
be redesigned. 

This method was applied to a nozzle design problem as shown in Figure 17. 
In that figure, the target geometry is shown as a solid line, and corresponds 
to the pressure distribution shown as a solid line in Figure 18. The 
pressure distribution was input to the design calculations as a target 
pressure, and the geometry shown as a broken line in Figure 17 was used as an 
initial guess. The best geometry, represented by the circles, was obtained 
after 15 wall corrections and 670 iterations of the direct solver, which 
shows that this method is somewhat slower than the previous method. The 
corresponding pressure distribution, shown by the circles in Figure 1 8 , is 
quite close to the target pressure. 

In order to check if results could be obtained with an initial geometry 
taken farther from the correct geometry than in the above case, the same 
problem was restarted with an initial geometry as shown in Figure 19. The 
best computed geometry, obtained after 22 wall corrections (700 direct solver 
iterations), is also shown in Figure 19. Here as before, the corresponding 
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Non-dimensional Axial location x/r 0 


Figure 17. Wall geometry at various iteration levels. Virtual velocity 
inverse method applied to a nozzle at Re<j=24. 



Non-dimensional Axial location x/r 0 


Figure 18. Wall pressure distribution prescribed versus computed. Virtual 
velocity inverse method applied to a nozzle at Re<j=24. 
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— Target geometry 



Figure 19. Wall geometry at various iteration levels. Virtual 

velocity inverse method applied to a nozzle at Re d =24 . 


Unmodified blade 



Figure 20. Wall mach number distribution at various iteration levels. 

Virtual velocity inverse method applied to a portion of a 
compressor cascade. Inlet mach number 0.44. 
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pressure was quite close to the target pressure. In the above cases, A = 0.2 
was used and the sign of the R.H.S. of equation 17 was taken as positive if 
p d was greater than p and negative otherwise. More illustration may be found 
m reference 30. 

3. 2. 2. 2. Redesign of Compressor Cascade Blade 

Since the relationship between pressure and wall geometry may become 
quite complicated in a cascade problem, the selection of two initial 
guess-geometries may not be as easy in this case as it is m nozzle or 
diffuser flows. For that reason, the virtual wall velocity method was 
selected for use in cascade computations. In order to test the method, the 
compressor cascade previously described for direct computations was again 
considered, with the inlet Mach number reduced to 0.44 so that transonic 
behavior could be avoided. The wall Mach number from a direct computation on 
the unmodified cascade is shown in Figure 20. For the design calculations, 
it was decided that the upper surface Mach number "bump", as seen in Figure 
20, be removed. The wall correction method was applied only in that portion 
of the blade, and the rest of the geometry was kept the same as for the 
unmodified blade. Such a modification only serves the purpose of 
illustrating the inverse design, and it is understood that in most cases, one 
may need to increase the area inside the Mach number curve rather than 
decrease it. 

For the cascade computations, A was not kept constant, but was 
calculated as 



where A' is specified. The desired effect of this modification was to 
emphasize the difference between small and large values of | Ar | calculated 
from Equation (17). Particularly, geometry corrections near the junctions 
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between the specified geometry region and the specified pressure region were 
to be made smaller, while corrections inside the specified pressure area 
(where | P d - P | is larger) were to be made comparatively larger. 

The design calculations were initiated using the unmodified cascade and 
its converged solution. After each wall correction, the new geometry was 
manually input to Sorenson's Grape code to obtain a new grid. It should be 
noted that making the process fully automatic requires coupling of the grid 
program with the Navier-Stokes code). After the 6th wall correction, X' was 
increased from 0.5 to 0.9 to accelerate the process. The final Mach number 
distribution shown in Figure 20 was obtained after 26 wall corrections (^,800 
iterations). It is quite close to the target distribution, and could have 
been made even closer by allowing for more wall corrections. The 
corresponding geometry modification is shown in Figure 21. 

Examination of the wall Mach number history after each wall correction 
reveals that the wall Mach number in the region being modified behaves in a 
wave-like manner. A wavy wall Mach number distribution such as the one shown 
in Figure 20 (10th wall correction) is typical. The waves are generated at 
the right side of the region being modified and move to the left, where they 
tend to vanish. The process will be stable if the amplitude of each new wave 
is smaller than the previous one. 

The speed of the calculations is controlled by both At' and X'. Since 
At' is usually limited by stability requirements for the direct solver, 
acceleration of the inverse process is obtained by varying X'. At this 
point, the upper limit for X' has not been established as yet. 

At the junctions between the region of specified geometry and the region 
of specified pressure, little control can be exercised over the pressure 
distribution, as shown by the Mach number spike on the left side of the 
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Figure 21. A portion of blade shapes before and after modification 


corrected region. Similar behavior has been reported in Reference 31. 

Removal of such a spike can be done by independent smoothing of the geometry 
or by extending the region of specified pressure. 

In the cascade computations, the sign of the right hand side of Equation 
(17) was taken as positive for P d < P and negative for P^ > P. The method 
works well for subsonic cases. However, for transonic calculations where 
small changes in geometry may lead to substantial changes in flow behavior, 
particularly when shocks are present, a more elaborate way of determining the 
sign of the wall correction may be required. 
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4. Conclusions 


f 

The following objectives of this research were achieved: 

1. An implicit factored algorithm was applied to the two-dimensional 
(axisymmetric and planar) thin layer Navier-Stokes equations. Good results 
were obtained for a variety of cases. 

2. Two methods for using the direct solver as an inverse design tool have 
been examined. The method of virtual velocity was used for compressor 
re-design calculations. In the re-design illustration, 26 wall corrections 
or 4800 time-steps were required to yield the final geometry. The rate of 
convergence is considered quite acceptable. 

5. Recommendations for Future Work 

It is believed that if one gains more experience in using the direct 
solver one may obtain more accurate results with less computation time. 

Also, improvements can be made for more efficiency and better accuracy. For 
example, adaptation of the algorithm to vectorized computers may lead to 
substantial savings in computer time. Introduction of upwind differencing 
would lead to a better handling of shock waves in transonic flows. For flows 
with supersonic inflow or outflow boundaries, conditions compatible with the 
supersonic nature of the flow will have to be imposed. Also, higher level 
turbulence models can be used, and extensions to the three-dimensional case 
are desirable. ' . 


o 
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APPENDIX A 


Body-Fitted Coordinate System 

In order to have boundary conditions satisfied exactly over an arbitrary 
shaped body, the Navier-Stokes equations are written in a transformed 
coordinate system such that the boundaries of the flow field correspond to 
grid lines of the finite difference network. In this Appendix, a method of 
obtaining such a transformation is explained. 

Following Thompson et al. (8), one may define a coordinate 
transformation from a physical plane (x, r) to a computational plane (£, n) 
so that the computational and physical coordinates are related by the Poisson 
equation: 

*«rr ‘ '<«•"> 

\x * n rr ’ 6(4 ' n) • <»' > 


In the above, f and g are forcing functions to be specified. In practical 
problems, it is often desired to define a computational grid work in a 
rectangular, uniform, transformed plane. Equations (Al ) are therefore 
rewritten so that £ and n become the independent variables: 


a x__ - 2 8x,_ + Y x + — „(fx_ + g x ) = 0 

££ £n nn £ 6 n 


a r - 2 8r + Y r + - (f r + g r ) = 0 

££ £n nn j 2 £ n 


(A2) 


where 


2 2 2 2 
a = x + r , Y - x_ + r r 
n n £ 5 


8 = x x + r_r , J = = C n - C n 

£ n £ n x_ r - x r r x r r x 

£ n n £ 
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and 


€ 


x 


= J r 

n 



\ ' ■ J r c • "r ■ Jx { <«> 

Equations (A2) are a set of elliptic partial differential equations. If 
a set of values for x and r are specified at the boundaries of the region of 
interest, these equations can be solved numerically by using, for example, an 
SOR (successive over relaxation) or an LSOR (line successive over relaxation) 
technique . 

As mentioned earlier, f and g must be specified. If one selects f = g = 
0, equations (A1) reduce to a set of Laplace equations, and the solution of 
Equations (A2) yields a highly uniform distribution of grid lines as viewed 
in the physical plane. If a non-uniform grid is desired, for example when 
one needs a fine grid near a wall, which gets coarser as one moves away from 
the wall, then one must find a method of specifying f and g. Steger and 
Sorenson (9, 11) using the exponential form of f and g suggested by Thompson, 

derived a method of specifying the forcing functions so that orthogonality 

\ 

and a specified spacing can be maintained near the boundaries. According to 
this method, f and g have the following form: 

fU.n) = f (£) e" an ♦ f_(0 e -c(r|max _n) 

1 ^ 


g(S,n) = g 1 (C) e bn + g 2 (0 e d(nmax b) ^ 


(AH) 


where f t , f 1( g lf g t are positive constants which must be calculated so that 
four geometric constraints are satisfied, and a, b, c and d are positive 
constants to be specified. 
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The third and fourth constraints are the same as the first two, except 

that they are applied to the boundary j“j m ax (outer boundary). They lead to 

two equations similar to (10) for x I. . and r 

" ' J - J max " 

Now considering the boundary j = 1 (n=0), it is clear, since r^x is 
large, that Equations (AH) reduce to 

f(e.o) = f^o 

g(C.O) = g^U) (All) 

Assuming that Equations (A2) apply at the boundary j=1 , one may obtain 
expressions for fi(£) and f 2 (£) by combining (A2) and (All): 

f, (« - [J(r n h, - x^hj)] j., 

g. U) = [J(-r.h 1 + x h )]. (A12) 

where 

h = [-J 2 (a x _ -26 x + Yx )] . 
i e,e, tin nn j = i 


b=J 


max 


h 2 - [-J (a -28 % 


In the above equations, x n and r n are obtained from Equations (A10).. 

A similar procedure applied to the outer boundary allows one to obtain 
expressions for f 2 (£) and g 2 (£). Then for specified values of a, b, c, d, 
f(5,n) and g(5.n) are known functions which can be used in the solution of 
Equations (A2). For details on the solution procedure, the rea'der is 
referred to Reference 9. 

In this research, a computer code was written to solve Equation (A2), 
both with f =g=0 and with f and g specified as explained above. The code was 
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tested on simple cases. However, for complicated viscous grids such as for 
airfoil cascade calculation Sorenson's Grape code is used (Ref. 9). 



APPENDIX B 


Matrices used In the Numerical Algorithm 
In this appendix, the matrices used in the algorithm for the thin layer 
Navier-Stokes equations are presented. Equations (12) and (13) of the text 
were derived for the numerical solution of Equation (1) of the text. All the 
matrices involved are listed as follows: 



0 

i 

-1 

1 

0 i 

0 

A = - 

3-Y 2 1 -Y 

2 U 2 

2 I 

v 1 

L_ 

(Y-3)u 


(Y-1)v 1 

1 

1-Y 


uv 

1 

| .. 

-V 

1 

1 

1 

0 


Yeu + (1_ Y ) 
P 

U (u2 + v2)j - 

+ II 

p 

" T' 

I(3u2+v2) 

2 1 

-p 

(Y-l)uv | 

-yu 


0 

1 

0 

1 _1 


0 " 


uv 

1 

I 

-v 

| 

1 “ U 


0 

B = - 

3H v 2 + HI 
2 2 

2 j 

u 1 

( Y 1 )u 

1 

| (Y-3)v 


1-Y 


— + (1-Y) 
P 

| 

v (u 2 + v 2 ) j 

( Y — 1 ) uv 

1 Ye^ Y- 

1 “ P 2 

i(3v 2 + a 2 ) 

-Yv 


From the definition of A and B we can show that 
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B is obtained from A by replacing by n x and by r^. 
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«3 - M (n x - T n r J 


a, = r yti n 
2 3 x r 


a Y |_ ( n 2 

4 Pr v x r ' 


a.r a.r a.r 

- + 2 f— ] + (—1—1 

[ J J i,j + 1 * L J J i,j 1 J J i,j-1 


a.r 


a.r 


‘1 


i.J + 1 




a.r 


a.r 


i J-1 


^ J ^ij + ^ J 


a 2 r 


a 2 r 


a 2 r 


1 2 iJ " ^ J \,j + 1 + J ^ i ,j + ^ J ^i,j-1 

a r a r 

1 2 iJ + 1 = HHi.J + 1 + HHi.J 

a r a r 

' (_:r)l ' J + (_r 


a_ is obtained by replacing the subscript 1 of a to 3 

3 i,J i,J 


is obtained by replacing the subscript 1 of a 


1 »j + 1 


i #j"1 


is obtained by replacing the subscript 1 of a 


i.J + 1 


i * j~1 


Similarly a^ , a^ and may be obtained. 


ij i,j + 1 


i.J-1 


to 3 
to 3 
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then 
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. 0 . . 
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0 

0 

0 1 

1 0 
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i »J + 1 
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i.j + 1 

1 ° 
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Matrices D, M 1 , E i and E i j_ 1 are given by 


D = 


0 0 

0 0 

o o 


0 

0 

-(Y-1 )v- 


4 u 


3 prRe 


0 

0 

Y-1 

0 


M. ■= — 
1 r 


2 v 


2 v 
■svn — 

3 rp 


4 v 

ttp — ( n u+n v) 
3 p x r 


0 

0 


2 v 


2 n x 


3 P 


2 n x 


3 P 

— (un +2vri ) 
3 p x r 


0 

0 

0 

0 


J i,j + 1 


±v_ 

3 Rej 




0 

0 


_0_ 

0 


0 

0 


, , ,Ju. 1 + 1 

\ i , j ( rp ) i ,j+1 I x ij rp 1,J 1 


+ ( n r ) i j Oi.J + 1 


{"VijWij+i 


Ij J _ 1 is obtained by changing i j + 1 of E,^ j + 1 to ^ j_ 1 . 
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APPENDIX C 


Solution of Block Tridlagonal Systems 
The use of the Beam and Warming implicit factored algorithm ( 23 ) in the 
solution of the thin layer Navier-Stokes equations leads to the inversion of 
tridiagonal systems. This appendix describes how the solution of such 
systems was obtained in this research. 

Consider the following system: 


Ql Ni 


~ X 1 “ 


~ F 1 " 

Lg Q2 ^2 


x 2 


F 2 

L3 Q3 N3 

• • • 

• • • 


x 3 

• 

• 


f 3 

• 

• 

• • • 

L n-1 Qn-1 N n-1 


• 

x n-1 


• 

F n-1 

L n Qn _ 


_ x n _ 


_ F n _ 


(Cl) 


as may result from Equations ( 12 ) or ( 13 ) of the main text. Here, L*, Qj , 
and Nj are 4 x 4 matrices, and and F ^ are 4 element column vectors. 

Equation (Cl) is rewritten as: 

EX = F , (Cla) 

where "=" denotes a matrix and a vector. A Lower-Upper decomposition of 
E yields: 
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where uj and lj are 4x4 matrices, and I is the 4x4 unity matrix. Carrying 
out the multiplication on the right hand side of Equation (C2) and 
identifying terms, one obtains: 


Hi = Q x 

(C3a) 

4 i = Li u^ 

(C3b) 

ui = Qi - J-i Ni_i 

(C3c) 
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with i = 2, 3, . . . n 

Now replacing E; in Equation (Cla) by its Lower-Upper decomposition, 
carrying out the multiplication and identifying terms, one obtains: 


gl = fi (C4a) 

and gi = fi “ 8>i gi-i , i = 2, 3, ... n (C4b) 

X n = u n ^ 8n (C5a) 

Xi = Ui -1 (gi - NiX i+1 ) , i = n-1 ... 2, 1 (C5b) 


Therefore, knowing Ui and 1^ from Equations (C3) and gi from Equations 
(C4), one may calculate the elements of X from Equations (C5). It should be 
noted that all required matrix inversions apply to 4x4 matrices and can be 
easily performed. 
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APPENDIX D 


Boundary Conditions Extrapolation 

The following sketch shows the boundary conditions specified for pipe 
flows . 


at j = j 


max 


A 3T 3p 

u = v = 0 -r- = -^ = 0 

3 n 3p 


p 

u 


specified 


J - jmax 


p specified 
j = 1 


i=1 


at j 


1 


v = 0 


3n 


_3T 

3n 


3p 

3n 


= o 


For the explicit part, the R.H.S. of Equation (10), at the inflow, the 
pressure p is obtained from linear extrapolation: 

P *.j = ° + p 2.J " Hi p 3J * 


AS2 


where Asi = s 2 -si , As 2 = 3 3~S2 and 3 is the length along a E, = constant line. 
The energy term, e, is obtained from the knowledge of p, u, v and p. At the 
outflow, u, v and p are obtained from an expression similar to the above 
extrapolation formula. At the wall, j = j max , the following backward second 
order approximation is used to obtain T and p at wall from interior points. 

A i ^1 »Jmax ^’Imax-I + ^ »Jmax-2 nfA s 2 

*1 - ^ * 0(4n) . 

1 «Jmax 


v. , 3T 

where, -t~- = 0 and -r— 
3n 3n 


0 for an orthogonal grid. 


At the centerline, or line of symmetry, v = u n = T^ = 0 are specified. Again 
= 0. The following second order approximation is used to obtain u, T and 
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P : 


Ml.. 


* ^ 1.2 
2Ari 


+ 0(An 2 ) 


For the implicit part of the algorithm. Equation (13) shows that at j = Jmax - ^* 

the knowledge of N . q? . or a typical term such as B^ 1 , Aci* 1 . is required, 
jm i,jm i.jm i » jm 

At the wall, u = v = 0 and the matrix multiplication of (i'Vq 11 ), . 

1 1 jm 

yields the following: 


(Y-1) 


0 

(n x Ae h )i,j m 
(nrAe^)i t jm 
0 


For p n ■= 0 one may approximate Ae^ t j m as follows. 


Ae i,jm = 3 ^i,jm tr Ae \ , j m ~ 1 3 ^A,jm A e A , jm-2 


r J — n 


1 rr' 




. n . — n 

= a Ae, , , + b Ae, , „ 

i , jm-1 i , jm-2 


then 


(B n Aq n ), . may be rewritten as 
i.jm 
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0 0 0 0 
0 0 0 (n x )i jm-i (Y-1)a 

0 0 0 (r »r>i 

0 0 0 0 




,jra-1 + 


0 0 0 0 
° 0 0 (r 'x)i,jm(Y-1)b ] 

° 0 0 ( ^)i, jm (Y-1)b! 

pooo 


Aq^ 


J 


In such a manner, N aTT 0 

Jm Aq i,jm 13 obta ined as a linear combination of 


i n _ n 

Q i Jm-i and Aq i , jm-2, Mhlch are the values at the interior 


points , 


,jm-2 


Th 



APPENDIX E 


Computer Program 

The computer program developed for the solution of the two dimensional 
(axl symmetric and planar) thin layer Navier-Stokes equations is presented in 
this appendix. The structure of the code is discussed, and a description of 
the code input variables is given. Output from the code is also described, 
and suggestions concerning how to run the code are given. 

Code Description 

The main routine of the code is very short, and merely allows to select 
between Equations (3), ( *0 or (5) of the text. It contains a parameter 
called IAXISY. For IAXISY = 0, Equation (4) is used. For IAXISY = 1, 
Equation (3) is used, and for IAXISY = 2, Equation (5) i3 used. For these 
three values of IAXISY, calls are made to Subroutines AXISYO, AXIS Y1 and 
AXISY2 , respectively. These three subroutines are similar to each other in 
their internal structure: they all have an initialization part, or solution 

part, a boundary condition part and a wall correction part. 

In the initialization part all constants and variables needed for the 
calculations are assigned values. A detailed description of input data is 
given below. During the initialization, a call is made to Subroutine TRANSF 
for the calculation of the finite difference grid. Also, during the 
initialization, initial values are given to the flow variables p, u, v, e, P 
(in the code, RH02 (I,J), U2(I,J), V2(I,J), E2(I,J), P2(I,J)) in the entire 
flow field considered. 

In the solution part, the right hand side of Equation (1) of the text 
(the Kj^j's) is formed, with the explicit artificial viscosity terms included 
The right hand side is stored in arrays HR(I,J), HM(I,J), HN(I,J) and HE(I,J) 
Then the £-sweep is carried out: matrices L^, Qj , and (represented by 
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C(I, K1, K2) , A( I , K1, K2) and B(I, K1, K2) , where K1 and K2 have values from 
1 to ty) are constructed, with the implicit artificial viscosity terms 
included, and the system of equations represented by Equation (12) of the 
text is solved. During the solution of Equation (12), a call is made to 
Subroutine TRIDIA which performs block tridiagonal inversions. Subroutine 
TRIDIA itself makes calls to Subroutine INVERT, which performs the inversion 
of !) x >1 matrices. The solution of Equation (12) yields the Aqfj * s, 
which are stored in HR(I,J), HM(I,J), HN(I,J) and HE(I,J) for use in the 
^-sweep. During the p-sweep, matrices Lj , Qj and Nj are constructed with 
implicit artificial viscosity terms included, and a call to Subroutine TRIDIA 
leads to the solution of Equation (13) of the text for the Aq^ t j ' s. 

In case the flow is turbulent, a call is made, at the beginning of the 
solution part, to Subroutine TURB or SUBROUTINE CEBECI, for the Baldwin-Lomax 
or the Cebeci turbulence model. 

In the boundary condition part, the values of the flow variables are 
obtained at the boundaries according to the guidelines given in the text. 

Only the explicit treatment is done in this part, since the implicit 
treatment is involved in the formation of matrices Lj, Q^, and Lj , Q j , Nj . 
In cases where scalar tridiagonal inversions are to be done (as in airfoil 
problems as explained in the text), Subroutine STRID is called. 

The wall correction part is only used when an inverse design is 
desired, and only if the wall pressure residuals become smaller than rp. 

When the secant method is used, a call is made to Subroutine SMOOTH, for the 
least squares smoothing of the computed geometries. 

Subroutine TRANSF is the grid generation subroutine. In this 
subroutine, the coordinates of the boundary points are input and the 
coordinates of the interior modal points (X(I,J), R(I,J)) are obtained in two 
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ways: by interpolation from boundary points or by solution of elliptic 

partial differential equations (Appendix A). In addition, the grid may be 
read from a separate file (i.e. when the grid has been calculated from a 
different code and stored in a file). It is also possible to obtain a 
coarse, but smooth grid by solving the equations of Appendix A and refine it, 
say near a solid wall, by a spline interpolation procedure. 

Input variables 

The following is a description of the input variables, their meanings, 
and how they should appear on data cards. A complete description is first 
given for the planar case. Then changes related to the axisymmetric case are 
shown. 

For the planar case, data are input as follows (M, L and all variables 
beginning with I or J are right justified integers): 


Card 1 

Column 

Variable 

Description 

1-10 

IAXISY 

This flag is zero for planar flow. It is one 
for axisymmetric flow if Eq. (3) is used, and 2 
if Eq. (5) is used. Eq. (3) should be used only 
for straight channels. Eq. (H) can be used for 
more complex shapes, but has not been tested for 
Re>1 00000. 

Card 2 

Column 

Variable 

Description 

1-10 

STEP 

Time step for the numerical algorithm. Input in 
the form X.XXXE+XX, right justified. 

11-20 

Z 

Numerical scheme selector. For the 
three- point-backward scheme, Z = 0.5. 

21-30 

THETA 

Numerical scheme selector. For the 
three-point-backward scheme, THETA = 1.0. 

31 - 40 

M 

Maximum number of grid points in the i-direction 
(^-direction) . M should be less than 80. 
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41-50 


L Maximum number of grid points in the j- direction 

(n - direction). L should be less than 45. 


Card 3 



Column 

Variable 

Description 

1-10 

UMAX 

Maximum number of iterations (time steps) 


allowed. 


11-20 

IFREQ 

Frequency with which intermediary data are 
output. If no intermediary output is desired, 
make IFREQ greater than UMAX. If output of 
the dependent variables in the entire flow field 
is desired every 10 iterations, IFREQ = 10. 

21-30 

CRIT1 

Value below which the density residual must fall 
for convergence. 

31-40 

CRIT2 

Value below which the residual of pu must fall 
for convergence. 

41-50 

CRIT3 

Value below which the residual of pv must fall 
for convergence. 

51-60 

CRIT4 

Value below which the residual of e must fall 
for convergence. 

61-70 

CRIT5 

Value below which the wall pressure residual 
must fall before a wall correction is made in 
design calculations. If no inverse design is 
desired, make CRIT5 smaller than CRIT1 , CRIT2, 
CRIT3, CRIT4. CRIT1 to CRIT5 are input in the 
form X.XXE-XX and right justified. 

Card 4 

Column 

Variable 

Description 

1-10 

P0 

Inlet stagnation pressure (Psia) 

11-20 

TO 

Inlet stagnation temperature (°R). 

21-30 

AMACH 

Inlet Mach number. For IFLAT=0, (see below), 
AMACH and the stagnation quantities are used to 
calculate the inlet velocity UINF, which is kept 
constant and is used as a reference velocity. 

The same thing applies when IFLAT=1 and IP0UT=0. 
For IFLAT=1 and IP0UT=1, this Mach number is 
discarded and recalculated from P0 and POUT. 

0 
•=r 

1 

ro 

RAD 

Characteristic length (i.e. chord length for 
airfoil, inlet radius for nozzle, diffuser, 
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41-50 

51-60 

61-70 

Card 5 
Column 
1-10 

11-20 

21-30 

31 _J *0 

41-50 

51-60 

61-70 

Card 6 
Column 
1-10 


etc ) in feet. RAD is input in the form 

X.XXXE+XX. 

ALPHA Artificial viscosity coefficient for the 

explicit part (e- 1 ). 

IV IS IVIS=1 for viscous flow. IVIS=0 indicates an 

inviscid flow. However, adjustments have not 
been made for inviscid calculations. 

IFLOW Use zero if only grid computations are needed, 

and one if both grid and flow computations are 
needed. 


Variable 

IREAD 

ITURB 

IFLAT 

INP 

IPRT 
IS YM 
ICASC 


Description 

Use zero if initial flow data are internally 
prescribed (usually uniform conditions) and one 
if initial flow data are read from a file. 

Use zero for laminar calculations and one for 
turbulent calculations. 

Use one for airfoil calculations, and zero for 
flow between parallel plates, symmetric nozzles 
or diffusers, etc., (where the centerline is a 
boundary) . 

For IFLAT=0 , use INP=0 if inlet conditions are 
internally prescribed (usually uniform 
conditions) and INP=1 if inlet velocity 
distribution is manually input. Not used if 
IFLAT = 1 . 

Use I P RT = 1 if initial flow data throughout the 
flow field are to be output, and IPRT=0 
ot her wise. 

Use ISYM=1 for symmetric airfoils at zero 
incidence angle and ISYM=0 otherwise. Not used 
for IFLAT=0. 

Use one for cascades, and zero for isolated 
airfoils. Not used for IFLAT=0. 


Variable Description 

IB Number of points along each periodic boundary in 

cascades. For isolated airfoils, IB=1 . 
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11-20 

IC 

IB + number of points along the inlet boundary 
+1. For isolated airfoils, IC = M. 

21-30 

IG 

Number of points along the wake boundary 
(between and including exit and trailing edge). 

31*40 

IGP 

IG + Number of points along the airfoil surface 
(counting the trailing edge only once). IB, IC, 
IG and IGP are disregarded for IFLAT = 0. 

Card 7 



Column 

Variable 

Description 

1-10 

ITEMP 

Use zero if temperature is not maintained fixed 


at the inlet, and one if it is (in which case 
velocity is also kept fixed at the inlet and 
pressure is extrapolated). Not used for 
IFLAT=0 . 

11-20 IPOTO Use zero if stagnation pressure and stagnation 

temperature are not kept fixed at the inlet, and 
one if they are (in which case flow angle is 
prescribed and Riemann invariant is 
extrapolated). To maintain p, u and v fixed at 
the inlet, use ITEMP=IP0T0=0. For IFLAT=0, use 
ITEMP=IP0T0=0. 

21-30 THIN Prescribed angle between the inlet flow 

direction and the x-axis (degrees). Not used 
for IFLAT=0. 

3 1 -40 THOUT Rough estimate of the exit flow angle (degrees). 

Used in internal calculation of initial velocity 
field. Not used for IFLAT=0. 

41-50 XTE The x~coordinate of the trailing edge, 

non-dimensionalized with chord length. Not used 
for IFLAT=0. 

Card 8 

Column Variable Description 

I- 10 ITWALL Use zero for an adiabatic wall and one if wall 

temperature is prescribed. Not used for 
IFLAT=0. 

II- 20 TWALL Prescribed wall temperature for ITWALL=1 (°R). 

21-30 POUT Prescribed exit pressure (Psia). Not used if 

IFLAT=1 and IP0UT=0. 
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31 -*10 


Card 3 
Column 
1-10 

11-20 

21-30 

31-40 

41-50 

51-60 

61-70 

Card 10 
Column 
1-10 

11-20 

21-30 


IPOUT 

Use IP0UT=1 if the exit pressure is input as 
POUT. Use IP0UT=0 if the exit pressure must be 

internally specified as the free stream 
pressure . 

Variable 

Description 

INVERS 

Use one if inverse calculation is desired, and 
zero otherwise. 

NEWAL 

For IFLAT=0, use NEWAL=0 for secant method and 
NEWAL=1 for virtual wall velocity method (for 
IFLAT=1 only virtual wall velocity method is 
used). 

IDES1 

Value of i corresponding to the beginning of the 
controlled region in inverse design. For 
IFLAT=0 use IDES1=1. 

IDES2 

Value of i corresponding to the end of the 
controlled region in inverse design. For 
IFLAT=0, use IDES2=M. For airfoil calculations , 
the controlled region must be on the upper 
surface, unless modifications to the code are 
made. For IFLAT=0 (and for the axisymmetric 
case), only the upper boundary (J=L) can be 
modified. 

FACT 

Multiplication factor used in smoothing for 
virtual wall velocity method (0.2 has been 
used) . 

RLOW 

Minimum allowable value of wall distance from 
centerline in wall correction calculations with 
secant method. Used in smoothing, to discard 
erratic values. 

RHIGH 

Maximum allowable value of wall distance from 
centerline (see RLOW). 

Variable 

Description 

IM 

Same as M for grid calculations where ICLUS1=0 
or when IS0REN=1 (see below). 

IM 

Same as L for grid calculations where ICALC=0 or 
ICLUST=0 or when IS0REN=1 (see below). 

UMAX 

Maximum number of iterations for grid 
calculations. 
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31-40 


I CALC 


41-50 

Card 11 
Column 
1-10 

Card 12 
Column 
1-10 


11-20 

21-30 


If ICALC=0, the grid is obtained by 
interpolation from specified boundary values. 
If ICALC=1 , a PDE solution (APPENDIX A) gives 
the grid. 

W Relaxation coefficient in SOR solution of grid 

PDE ' s ( 0< VK 2 . 0 ) . 


Variable Description 

IS0REN If IS0REN=1 , the grid is read from a file (e.g. 

a file obtained from Sorenson's program). If 
IS0REN=0, the grid is internally calculated. 

(skip if IS0REN=1 ) 


Variable Description 

ICLUST For simple channels (diffusers, nozzles, etc...) 

the grid can be obtained as follows: calculate 

a coarse but smooth grid (ICALC=1) and then 
refine the grid by adding grid lines, say near 
the wall. If ICLUST-0, no grid refinement is 
needed. If ICLUST=1 , the coarse grid is 
calculated using the input values of X(1,J) and 
R(1,J) (inlet coordinates). Then for the finer 
grid, a new distribution of the X(1,J)'s and 
R(1,J)'s must be given (the XG(J) and RG(J) 
values). With the new inlet distribution, the 
maximum number of points in the J-direction 
changes from JM to JMC. The new distribution of 
inlet points does not have to cover the entire 
inlet. For example, one may need to keep points 
from J=1 to J=JMLL as they are, and input a ne,w 
distribution for J=JMLL+1 to J=JMC (see Figure 
E-1). The XG' s and RG's are used to read in 
these JMC-JMLL values, and the grid in the rest 
of the domain is obtained by spline 
interpolation from the old grid. 


JMC Maximum number of points in J-direction for 

refined grid when ICLUST=1 is used. For 
ICLUST=1, JMC is equal to L. 

JMLL Number of inlet grid points (starting with J=1) 

which are left unchanged when a refined inlet 
distribution of points (ICLUST=1) is prescribed. 
If the refined distribution covers the whole 
inlet, make JMLL=1 . 
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• Coarse distribution (X(I,J),R(I,J)) 



Figure E-1 : Grid refinement option (ICLUST=1) 

i 


Card 13 (Skip If IS0REN=1) 


Column 

Variable 

Description 

1-10 

ICLUS1 

Similar to ICLUST, but applies to I-direction. 
Associated with IMC, IMLL, XGI(I), RGI(I). The 
ICLUS1 feature is not used for inverse design 
computations where grid computations are 
performed many times during the same run, after 
each wall correction. 

11-20 

IMC 

Similar to JMC, but applies to I-direction. For 
ICLUS1=1, IMC is equal to M. 

21-30 

IMLL 

Similar to JMLL, but applies to I-direction. 

Card 1 4 

(Skip If IS0REN=1 ) 


Column 

Variable 

Description 

1-10 

PI 

Initial value of f-| (Eq. A4). A value of zero 
is usually used. 

11-20 

PX1 

Coefficient a in Eq. (A4). 

o 

on 

c\j 

P2T 

Initial value of ( E Q* Al< )* A value of zero 

is usually used. 

31-40 

PX2 

Coefficient c in Eq. (A4). 
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Card 15 (Skip if IS0REN=1) 


Column 

Variable 

Description 



1-10 

Q1 

Initial value of g-j in Eq. (A4). A 
zero is usually used. 

value 

of 

11-20 

QX1 

Coefficient b in Eq. (A4). 



21-30 

Q2 

Initial value of gg in Eq . ( A4) . A 
zero is usually used. 

value 

of 

31-40 

QX2 

Coefficient d in Eq. (A4). 



Card 16 

(Skip if IS0REN=1 ) 




Column 

Variable 

Description 



1-10 

WP 

Relaxation parameter in solution of 
equations (0<WP<2.0) . 

grid 


11-20 

WQ 

Similar to WP 



21-30 

PLIM 

Limitation factor to help stability of grid 
solution (PLIM=1.0 has been used). 

o 

m 

QLIM 

Similar to PLIM. 



Card 16 

(Skip if IS0REN=1 ) 




Column 

Variable 

Description 



1-10 

WR 

Similar to WP 



11-20 

WS 

Similar to WP. 



21-30 

RLIM 

Similar to PLIM 



31-40 

SLIM 

Similar to PLIM. 



Next input 





. Values of X(I f 1) f the x values of grid points along the boundary J=1 . IM 
values should be input, 5 per card, in 13 columns fields. Data are input in 
the form X . XXXXXXE+XX . All values are non-dimensionalized by characteristic 
length and are right justified. Skip if IS0REN=1 . 
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Values of R(I,1), the r values of grid points along the boundary J=1 . 

IM values must be input. Skip if IS0REN=1 . 

Value of X( I , JM) , the x values of grid points along the boundary J=JM. 

IM values must be input. Skip if IS0REN=1 . 

Values of R(I,JM), the r values of grid points along the boundary J=JM. 
IM values must be input. Skip if IS0REN=1 . 

Values of X( 1 , J) , the x values of grid points along the 1=1 boundary. 
JM-2 values are needed (J=2 to J=JM-1). Skip if IS0REN=1 . 

Values of R ( 1 , J ) the r values of grid points along' the 1=1 boundary. 
JM-2 values are needed (J=2 to J=JM-1). 

Values of X(IM,J), the x values of grid points along the I=IM boundary. 
JM-2 values are needed (J=2 to J=JM-1). Skip if IS0REN=1 . 

Values of R(IM, J) , the r values of grid points along the I=IM boundary. 
JM-2 values are needed (J=2 to J=JM-1). These values are used for the 
initialization but are not kept during the calculations. Skip if 
IS0REN=1 . 

Values of XG(J), the x values of nodal points along the inlet boundary 
for refined grid. JMC-JMLL values are needed. Skip if ICLUST=0 or if 
IS0REN=1 , or ICALC=0. 

Values of RG(J), the r values of nodal points along the inlet boundary 
for refined grid. JMC-JMLL values are needed. Skip if ICLUST=0 or if 
IS0REN=1 , or ICALC=0. 

Values of XGI(I), the x values of nodal points along the J=1 boundary 
for refined grid. IMC-IMLL values are needed. Skip if ICLUS1=0 or 
IS0REN=1 or ICALC=0. 
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Values of RGI(I), the r values of nodal points along the J=1 boundary 
for refined grid. IMC-IMLL values are needed. Skip if ICLUS1 = 0 or 
IS0REN=1 or ICALC=0. 

Values of PDES(I), the target pressure distribution, in Psia, for 
I=IDES1 to I=IDES2 (from 1=1 to I=M for axisymmetric case). Skip if no 
design calculations are desired. 

Values of U2(1,J), the inlet velocity distribution, for J=1 to J=L, in 
ft/sec. Skip if INP=0. 

Values of the r coordinate at the wall (R(I,JM)) for second guessed 
geometry if secant method is used. Data are input for 1=1 to I=M. 

Skip if secant method is not used. 

In the above, all boundary data are input 5 per card, in the form 
x.xxxxxE+xx, and are right justified. Geometry data are 
non-dimensional. 

For the axisymmetric case, data are input in the same way, except 
for the following cards: 


Card 3 



Column 

1-10 

Variable 

PINF 

Description 

Exit static pressure, used as initial pressure 
everywhere (Psia) 

11-20 

TINF 

Initial temperature, used as reference 
temperature (in °R). 

21-30 

UINF 

Inlet velocity, used as reference velocity (in 
ft/sec) . 

31 -40 

RAD 

See planar case. 

Ml -50 

ALPHA 

See planar case. 

51 -60 

IVIS 

See planar case 

61-70 

IFL0 

Same as IFL0W in planar case. 
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Card 4 


Column 

Variable 

Description 

1-10 

ITURB 

See planar case. 

11-20 

INP 

See planar case. 

21-30 

IREAD 

See planar case 

0 

■=T 

1 

cn 

IPRT 

See planar case. 

Card 5 



Column 

Variable 

Description 

1-10 

ITEMP 

Use ITEMP=0 if density is kept fixed at the 
inlet, and ITEMP=1 if temperature is kept 
constant. In both cases the inlet velocity is 
kept fixed and pressure is obtained by 
extrapolation. 

11-20 

TO 

This stagnation temperature and UINF are used t< 
obtain the static temperature which is 
maintained fixed at the inlet. 

21-30 

ICENT 

Use ICENT=1 if the centerline is a boundary of 
the flow field (J=1) and ICENT=0 if there is a 
solid wall at J=1 (annular type). For ICENT=1, 
only laminar calculations can be performed, 
unless subroutine CEBECI is modified to account 
for two walls. 

Card 6 



Column 

Variable 

Description 

1-10 

INVERS 

See planar case 

11-20 

NEWAL 

See planar case 

21-30 

FACT 

See planar case 

31-110 

RL0W 

See planar case 

41-50 

RHIGH 

See planar case 


For the axisymmetric case, skip cards 7 to 9 of the planar case. 
Output Description 

The output of the code first shows all the input parameters in the 
order in which they are input. Reference quantities such as velocity, 
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density, temperature, viscosity, etc which are calculated by the code, 

are also printed. 

Once the input parameters have been printed, the coordinates of the 
nodal points in the finite difference network are output: the X(I,J)'s are 

printed first, as M sets of data, each set corresponding to values of J from 
J=1 to J=JM. Then the R(I,J)’s are printed in the same manner. If ICALOO, 
these values are the ones obtained by interpolation. If ICALC=1 , these are 
the initial data for the grid computations. For ICALC=1 , the residuals DELX 
and DELR for X and R are then printed for the first and the last iterations 
of the grid computations. Then the final distributions of the X(I,J)'s and 
the R(I,J)'s are printed. 

If wall design calculations are to be done, the non-dimensional design 
pressure distribution (along the controlled portion of the solid wall to be 
modified) is printed next. If no design calculations are to be performed, 
the above is skipped and the initial velocity, density, pressure and energy 
fields are printed next if IPRT=1 . The form of the output is the same as the 
form of the X(I,J) and R(I,J) output. The printed values are 
non-dimensional, and the title before each set of data indicates how the 
non-dimensional values can be used to obtain the dimensional values. 

The residuals DELIP, DEL2P, DEL3P, DEL4P and DEL5P are printed next. 
Residuals DELIP to DEL4P are defined as the maximum differences between the 
values of p, pu, pv or e, respectively, between the current and the previous 
iteration in time (ITIME). DEL5P is defined the same way, but only applies 
to the wall pressure. Also printed are the locations (values of I and J 
such as I HR and JHR for p, IHM and JHM for pu, IHN and JHN for pv, IHE and 
JHE for e) where the maximum difference occurs. The residuals are printed 
for each iteration, and their overall trend should be decreasing. 
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If wall corrections are performed, the convergence history given by 
the above residuals will be interrupted whenever DEL5P becomes smaller than 
CRIT5. At this point, the wall pressure in the controlled region 
(corresponding to the current geometry) is printed. A new wall geometry is 
obtained, and the corresponding r-coordinates are printed. Also printed is 
DEL6P, the maximum difference between the current pressure and the target 
pressure. DEL6P is printed after every wall correction. However, it 
usually does not give an adequate indication on whether or not the current 
pressure is approaching the target pressure, since the wall correction scheme 
may not perform adequately at some points (such as the end points of the 
controlled region). A plot of the pressure distribution in the controlled 
region after each wall correction would give a better indication on the 
success of the wall correction scheme. For airfoil calculations, the 
calculations are stopped after each wall correction because a new grid must 
be externally obtained. 

When DELIP, DEL2P, DEL3P and DEH1P become smaller than CRIT1 , CRIT2, 
CRIT3 and CRIT 1 ), respectively (or after the number of iterations ITIME 
becomes larger than UMAX), the calculations are stopped, and the final 
density, velocity, pressure and energy fields are output. 

In addition to the output described above some data are automatically 
stored by the code on files. At the end of the grid calculations, the grid 
is stored in a data set with the reference number 13. At the end of the flow 
calculations, the entire flow field values of u, v, p, e, pu, pv, P are 
3tored in a data set with the reference number 12. This way, the grid and 
flow data can be read from data sets 13 and 12 and used, for example, in a 
plotting routine written by the user. 
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Additional Considerations 


The present code has been run in FORTRAN IV on an IBM 3081 computer. 
Double precision is used in the calculations. Adequate job control language 
statements must be provided in order to run the code. 

As mentioned in the input data description section, a grid file is 
read by the code if IS0REN=1 and a flow field data file is read if IREAD=1 . 
These files correspond to data set reference numbers 8 and 11, respectively. 
Adequate job control language statements must be incorporated to account for 
these input data sets. Also, in cases where the data sets are not available 
(which may be the case when IREAD=0 or IS0REN=0), if may still be necessary 
to create them (even if they contain no useful information) in order to make 
the execution possible. 

As mentioned in the output description section, some data are stored 
In data sets 12 and 13. Adequate job control language statements must be 
incorporated to account for these output data set3. One interesting way of 
using data sets 11 and 12 would be, for example, to start the computations 
with IREAD=0. The final values of u, v, p, e, pu, pv for the calculation are 
automatically stored in a data set named DATA1 for example, corresponding to 
reference number 12. If more iterations are needed for convergence, one may 
restart the calculations using IREAD=1 and, through the appropriate job 
control language statements, associate data set DATA1 to reference number 11 
and use a new data set name, i.e. DATA2, in association with reference number 
12. Then for the new run, the initial flow field is read from 11 and the new 
results are stored in 12 (DATA2). This allows for a more efficient use of 
the code in very long calculations. 

In order to read flow data stored by the code in a data set, a 
statement such as the following must be used: 


90 


c-a- 



DO SN1 J=1 ,L 
DO SN1 1=1 ,M 

SN1 READ(11 ,SN2), U2(I,J), V2(I,J), RH02(I,J), E2(I,J), M2(I,J), 

1 N2( I , J) , P2( I, J) 

SN2 FORMAT (5X, 7E1 5.6) 

In the above, SN1 and SN2 are statement numbers . The same format must be 
used in writing data into a file to be used by the code if IREAD=1. In order 
to read grid data stored by the code in a data set, a statement such as the 
following must be used: 

DO SN1 J=1 , L 
DO SN1 1=1 ,M 

SN1 READ(8,SN2) , X(I,J), R(I,J) 

SN2 FORMAT ( 2E1 4.6) 

Sample Input and Output 

The following is a simple example of a set of input data, and a 
computer plot of the velocity field (see Figure E~2). Notice that the plot 
was obtained by a small, separate code. Since plotting software varies with 
facilities, no plotting routine was incorporated in the code. 


« 
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2 

0. 500E+00 0.5 1.0 17 19 

450 1001 0.10E-04 0.10E-04 0.10E-04 0.10E-04 0.10E-05 

14.70 540.0 500.0 0.335E-02 2.0 1 1 

0 0 0 1 

1 542.0 1 

0 0 0.0 0.0 0.0 

17 11 100 1 1.55 

0 

1 19 1 

0 0 1 


0.0 

0.3 

0.0 

0.7 

0.0 

0.3 

0.0 

0.7 

0.3 

0.3 

1 .0 

1.0 

0.3 

0.3 

1 .0 

1.0 


0 . 00000E+00 0.12500E+00 0.25000E+00 0.37500E+00 0.50000E+00 
0.62500E+00 0. 75000E+00 0.87500E+00 0. 10000E+01 0.11250E +01 
0 . 1 2500E+01 0 . 1 3750E+01 0.15000E+01 0.16250E+01 0.17500E+01 
0 . 1 8750E+0 1 0 . 20000E+0 1 

0 . 00000E+00 0 . OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO 0.00000E+00 
0 . 00000E+00 O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO 
O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO 
O.OOOOOE+OO O.OOOOOE+OO 

O.OOOOOE+OO 0.12500E+00 0.25000E+00 0.37500E+00 0.50000E+00 
0 . 62500E+00 0.75000E+00 0.87500E+00 0.10000E+01 0.11250E+01 
0. 12500E+01 0. 1 3750E+01 0.15000E+01 0.16250E+01 0.17500E+01 
0. 18750E+01 0. 20000E+01 

0. 10000E+01 0. 10000E+01 0.10000E+01 0.10079E+01 0.10210E+01 
0. 10342E+01 0. 10473E+01 0.10604E+01 0.10736E+01 0.10867E+01 
0. 10999E+01 0.111 30E+01 0.11261E+01 0.11393E+01 0.11472E+01 
0.11 472E+01 0.11 472E+01 

O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO 
O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO 
0 . 1 OOOOE+OO 0 . 20000E+00 0.30000E+00 0.40000E+00 0.50000E+00 
0 . 60000E+00 0 . 70000E+00 0.80000E+00 0.88000E+00 0.95000E+00 
0. 20000E+01 0 . 20000E+01 0.20000E+01 0.20000E+01 0.20000E+01 
0 . 20000E+01 0 . 20000E+0 1 0.20000E+01 0.20000E+01 0.20000E+01 
0 . 1 OOOOE+OO 0 . 20000E+00 0.30000E+00 0.40000E+00 0.50000E+00 
0 . 60000E+00 0 . 70000E+00 0.80000E+00 0.88000E+00 0.95000E+00 
O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO 
O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO 
O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO 
O.OOOOOE+OO O.OOOOOE+OO O.OOOOOE+OO 

0 . 1 5473E+00 0 . 30946E+00 0.45157E+00 0.56525E+00 0.65620E+00 
0 . 72896E+00 0.78717E+00 0.83374E+00 0.87099E+00 0.90079E+00 
0. 92463E+00 0.94371E+00 0.95896E+00 0.97117E+00 0.98094E+00 
0.98875E+00 0.99500E+00 0.10000E+01 
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